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Abstract 

A strong-coupling expansion is applied to the anharmonic Holstein model and 
to the Holstein-Hubbard model through fourth order in the hopping matrix 
element. Mean-field theory is then employed to determine transition temper- 
atures of the effective (pseudospin) Hamiltonian. We find that anharmonic 
effects are not easily mimicked by an on-site Coulomb repulsion, and that an- 
harmonicity strongly favors superconductivity relative to charge-density-wave 
order. Surprisingly, the phase diagram is strongly modified by relatively small 
values of the anharmonicity. 
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I. INTRODUCTION 



The interaction of the conduction electrons in a solid with the lattice vibrations is de- 
scribed by the so-called electron-phonon problem. Mig dali and Eliashbergi pioneered the 
study of such interacting fermion-boson systems in the limit where the phonons are all har- 
monic, and the phonon energy scale is much smaller than the electronic energy scale. Vertex 
corrections can be neglected in this case0 and a self-consistent theory can be constructed 
that is exact in the limit of weak-coupling; the theory is an expansion in powers of the 
coupling strength multiplied by the ratio of the phonon energy scale to the electronic energy 
scale. 

But phonons in real materials are never purely harmonic — higher-order (anharmonic) 
contributions to the phonon potential are always present. This phonon anharmonicity is 
responsible for many different physical effects in solids. For example, thermal expansion 
arises purely from anharmonic effects — a harmonic crystal does not change its volume upon 
heating. Such anharmonic effects have been treated in an approximate fashion: in the 
quasiharmonic approximation!! only the effect of thermal expansion is taken into account by 
postulating that the phonon frequencies have a dependence upon the volume of the crystal, 
as described by the Griineisen parameter; in the self-consistent harmonic approximation^ 
the harmonic force constants are self-consistently replaced by their thermal averages over 
all possible motions of the other atoms — it is used to describe systems with strong anhar- 
monicity; finally, in the pseudoharmonic approximation! both thermal expansion effects and 
phonon-phonon interactions are taken into account by employing the quasiharmonic approx- 
imation plus a perturbative expansion in the phonon-phonon interaction. Superconducting 
transition temperatures have also been studied, with the result that anharmonicity does 
not enhance the transition temperature in the weak to moderate coupling regime where 
Eliashberg theory appliesfl. But an exact treatment of lattice anharmonicity is difficult from 
a theoretical point of view because an anharmonic "perturbation" is never a small pertur- 
bation; the phonon wavefunctions are always dominated by the anharmonic terms in the 
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potential as the phonon coordinate becomes large. Even the qualitative effects of lattice 
anharmonicity on superconductivity are not well understood. All that is known rigorously 
about the anharmonic electron-phonon problem is that the ground state must contain a spin 
sinrietB for even numbers of electrons on a finite lattice. 

Much progress can be made however, in the limit of strong coupling (the electron-phonon 
interaction is much larger than the hopping integral of the electrons), where the electrons 
strongly bind together into preformed pairs (called bipolarons) and the ground state of the 
system is highly degenerate. Degenerate perturbation theory (in the kinetic energy of the 
electrons) about this bipolaronic ground state produces an expansion in inverse powers of 
the coupling strength. This theory has been exhaustively analyzed to second order in the 
hopping! and has recently been studied to fourth order in the hopping!. In this contribution, 
the fourth-order calculations are extended to include both anharmonic phonons and a direct 
electron-electron repulsion. 

The simplest electron-phonon model that includes both anharmonic effects and direct 
electron-electron repulsion is the anharmonic Holstein- Hubbard model00 in which the con- 
duction electrons interact with themselves and with local phonon modes: 



where c ia [c ia ) creates (destroys) an electron at site i with spin a, n ia = c ia c i(J is the electron 
number operator, and X{ (pj) is the phonon coordinate (momentum) at site i. The hopping 
of electrons between lattice sites % and j is governed by the hopping matrix element (t^ 
is a Hermitian matrix). 

The local phonon has a mass M and a frequency Q associated with it; the combination 
n := MVL 2 is a spring constant that measures the stored energy per unit length squared 
in the phonon coordinate. The anharmonic contribution to the phonon potential energy is 
chosen to be a quartic in the phonon coordinate with a strength a an h- The electron-phonon 
interaction strength is parameterized by an energy per unit length and is denoted g. A useful 
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combination of fundamental parameters is the bipolaron binding energy (in the harmonic 
limit with U c = 0) 

U--^—-- 9 - (2) 

which determines the energy scale for the effective electron-electron interaction mediated by 
the phonons. This attraction competes with the direct Coulomb repulsion denoted by U c . 
The chemical potential is \i. 

The hopping matrix elements t%j are used to define the energy scale. The choice tuvAr ='■ 
t* /2\fd with t* — 1 and d the dimensionality of the lattice, is made so as to have a well- 
defined limit when d — > oo. The mass is then set equal to one (M = 1) leaving U, U c , a an h, 
and Q as free parameters. The strong-coupling expansion is a perturbative expansion in the 
hopping terms of Eq. (JJ) and is valid when the bipolaron binding energy is much larger than 
the electronic hopping integral. Another important parameter in the strong-coupling limit 
is the polaron band-narrowing parameter denoted by S = \U\/Q. 

The original Holstein Hamiltonian corresponds to the case a an h = U c = 0. Both the har- 
monic Holstein model and the harmonic Holstein-Hubbard model have been solved exactly 
in the limit of infinite dimensions via quantum Monte Carlo simulationlllll!. These models 
display charge-density-wave (CDW) order near half filling and superconductivity (SC) away 
from half-filling. As the phonon frequency is increased, the SC is favored relative to the 
CDW order. However, in the strong-coupling limit, CDW order is favored over SC because 
of the band-narrowing effect of the bipolaron. The quantum Monte Carlo simulations also 
found that the effective phonon potential (determined after integrating out the effects of the 
electrons) generically acquires a double-well structure signifying the formation of a bipo- 
laron, and indicating that a strong-coupling expansion should be accurate even down to 
moderate values of the coupling strength. This has proven to be true for the harmonic casein 
and is likely to also hold over a more restricted region for the anharmonic case. 

An initial analysis of the anharmonic model in the strong-coupling limit (t^ = 0) can be 
made by using the Born-Oppenheimer approximation^: the phonon frequency is assumed to 



be smaller than any of the other energy scales so the phonons can be approximated by static 
lattice distortions corresponding to the minimum of the phonon potential energy. Since 
the phonons couple linearly to the electronic charge, the equilibrium phonon coordinate 
varies when there are zero, one, or two electrons on a site. The origin xq = is chosen to 
correspond to the case with no electrons on a site. Then x\ and £2 denote the equilibrium 
coordinates with one or two electrons on a lattice site. In the harmonic case the relative 
distances x\ — xq and X2 — x\ are identical, which is a requirement for particle- hole symmetry. 
When a lattice anharmonicity is turned on, the equilibrium phonon coordinates with one 
and two electrons on a lattice site all move toward the origin, but the relative distances are 
no longer symmetric, rather the distance x<i —x\ becomes significantly smaller than x\ — Xq, 
as can be seen by calculating the perturbative shift in the equilibrium coordinates^ as a 
function of the anharmonicity a an h 

_ n _ 9 4:a anh g 3 _ 2g | 32a anh g 3 

X °~ ' Xl ~ Mft 2 + AW ' X2 ~ Mft 2 + ilW • 1 ' 

This asymmetric shift of the equilibrium phonon coordinate as a function of lattice 
anharmonicity causes two main effects: (1) the model loses particle-hole symmetry which 
allows a new type of superconductivity to emerged and (2) the effective electron-electron 
attraction is sharply reduced as can be seen by a plot of the bipolaron binding energy 
in Figure 2 (b). Thus the lattice anharmonicity generates an effective retarded repulsive 
interaction between the electrons and breaks particle-hole symmetry, removing the nesting 
instability of the CDW at half filling and weak coupling. One expects that the lattice 
anharmonicity thereby to favor SC relative to CDW order, although it is also likely that 
anharmonic effects will reduce the transition temperatures (except close to the filled band, 
where the new hole-superconductivity mechanism can take over). What is surprising is that 
a rather small lattice anharmonicity can have a large effect on the electron-phonon problem. 

Since the main effects of lattice anharmonicities are driven by the asymmetric distri- 
bution of the equilibrium phonon coordinate when there is zero, one, or two electrons on 
a lattice site, one expects that anharmonic effects will be strongest in the small-phonon- 



frequency limit. This is because the phonon coordinates all approach zero in the high- 
frequency limit (because the phonon reacts instantaneously to the change in the electrons) 
and these asymmetric effects disappear. Since phonon frequencies tend to be small in real 
materials, anharmonic effects can be important even if the phonon potential energy appears 
to be well approximated by a harmonic potential, i. e., if a an h is small. 

Hirsch's new mechanism for superconductivity arises from an examination of the an- 
harmonic electron-phonon model in the static limitS The tunneling matrix element for a 
polaron from one lattice site to its nearest neighbor depends exponentially on the difference 
x\ — xq if there is no electron on the neighboring site and exponentially on X2 — x± if there 
is an electron occupying the neighboring site. The exponential dependence arises from the 
Franck-Condon overlap factors. Since these two values can be significantly different in an an- 
harmonic model, one finds that the electronic motion will be dominated by a kinetic energy 
that depends on the density of the electrons at a given site. This is precisely the physical 
situation needed to generate superconductivity from kinetic energy effects — paired electrons 
have a lower effective mass than a single electron, which induces the superconductivity tran- 
sition at a low enough temperature. This novel mechanism for generating superconductivity 



shown to exist in the exact solution of any model system. This new mechanism for super- 
conductivity disappears, however, in the strong-coupling limit where restriction is made to 
consider only empty sites and bipolarons. The single-electron (polaronic) states are inte- 
grated out because the bipolaron binding energy is much larger than the electronic hopping 
integral. Electron-hole symmetry is restored since these polaronic states only appear in 
virtual processes. 

The anharmonicity has a much different effect than a direct electron-electron repulsion. 
The electron-electron repulsion uniformly reduces the bipolaron binding energy without 
changing the equilibrium phonon coordinates when there are zero, one, or two electrons 
at a lattice site. Thus, (1) the bipolaron binding energy can become negative (signifying 
there is no electron-electron pairing) and (2) the system explicitly retains it's electron-hole 



has been analyzed via weak-coupling mean-field theory 
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symmetry. Furthermore, since the repulsion is instantaneous, the retardation effects are 
unchanged from the case without Coulomb repulsion, i. e., the Franck-Condon overlaps 
remain the same. Thus we expect CDW ordering to always survive at half filling (if there 
is a net electron-electron attraction), and that both SC and CDW order will disappear as 
the Coulomb repulsion becomes too large, but it is not clear whether or not they disappear 
together, or at different values of the electron-phonon coupling. 

This contribution is a continuation of the work of one of the authors! on the strong- 
coupling expansion for the harmonic Holstein model to include both anharmonic and 
Coulomb repulsion effects. In Sec. II, the formalism for the perturbation theory, and 
the generation of the effective pseudospin Hamiltonian will be described. In Sec. Ill, a 
mean-field-theory analysis of the pseudospin Hamiltonian will be given, and appropriate 
phase diagrams calculated for both the anharmonic case and the Coulomb repulsion case. 
Conclusions and a discussion will follow in Sec. IV. 



II. FORMALISM FOR THE PERTURB ATIVE ANALYSIS 

The strong-coupling expansions are carried out with a method based on perturbation the- 
ory. The ground state is a bipolaronic state consisting of either paired electrons or empty 
sites. The distribution of these bipolarons is not determined to zeroth order in the electronic 
kinetic energy, so the ground state is highly degenerate. The effective Hamiltonian (within 
this degenerate subspace) can be determined by using operator methods. In the late 1950 's, 
Andersonlll first used such methods to show that the strong-interaction-strength limit of 
the Hubbard model is described by a Heisenberg antiferromagnet with an exchange integral 
j = 4\Uj\ 2 /\U\ that vanishes as the interaction strength increases. Kato0 described how to 
determine the effective Hamiltonian for an arbitrary degenerate subspace using perturbation 
theory and operator methods. His analysis was applied to the Hubbard model in one dimen- 
sion by Klein and SeitzEl and in arbitrary dimensions by TakahashiEl (see also the recent 
work by van Dongen0). Beni, Pincus, and KanamoriS and Hirsch and FradkinS applied 
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the same methods to the harmonic Holstein model determining the effective Hamiltonian to 
second order in the hopping. This analysis was extended to fourth order by Freericksll. 

Kato's method^ begins with a Hamiltonian H = Hq + T with Hq the unperturbed 
Hamiltonian and T the perturbation. In our case, Hq corresponds to the Hamiltonian in 
Eq. (|I|) with t{j = 0, and the perturbation T is the electronic kinetic energy. The ground- 
state energy is Eq, Qq denotes the subspace that contains all of the degenerate ground states, 
and the projection operator onto Qq is Po' 

H P = P H = E P , P 2 = P • (4) 

As the perturbation is turned on, the eigenstates will evolve into a new subspace Q with 
corresponding projection operator P. If it is assumed that the subspace Q has a nonzero 
overlap with the unperturbed subspace Qq, then the standard eigenvalue equation (H — 
E)\E) = can be projected onto the unperturbed subspace Q Q : Pq(H — E)PP Q \E) = 0, 
to yield an effective equation for the perturbed eigenvalue E. The Hamiltonian PqHPPq 
acts purely within the unperturbed subspace Qq and has an overlap operator PqPPq that is 
not equal to the identity. Taking into account this nontrivial overlap, results in an effective 
Hamiltonian of the formll'i H e ff = Hq + H2 + H4 + . . ., since only even powers of the 
perturbation enter for the generalized Holstein model. The first two nontrivial terms of the 
effective Hamiltonian satisfy 

H 2 :=P T-^-TP , (5) 



and 



Ha . = p^ T 1 ~ P rp 1 ~ P rp 1 ~ P T p^ 

Eq — Ho Po — Ho Eq ~ H 



1 p 1 p 1 p 1 p 

P T- °-^TP T -TP + P T —T- \^TP 

(E — H ) Eq — Hq Eq — Hq (Eq — Hq) 



(6) 



The expansion for the effective Hamiltonian can be expressed graphically by a set of 
diagrams. A solid line denotes virtual processes where an electron hops from site i to site 
j with strength t^. All diagrams must be closed, since the effective Hamiltonian acts solely 
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within the degenerate subspace Qo implying each virtually broken electron pair must be 
restored. There is only one possibility for the second-order term, which corresponds to either 
hopping from site i to site j and hopping back to site i, or which corresponds to subsequent 
hops from site % to site j. The diagram that illustrates both of these processes is depicted in 
Figure 1(a). There are four possible diagrams that contribute to fourth order which are also 
depicted in Figure 1. The first three diagrams are linked diagrams which form nonvanishing 
contributions to the effective Hamiltonian. The last diagram 1(e) is an unlinked diagram 
which does not contribute to the effective Hamiltonian because the contributions from the 
positive and negative terms in Eq. (|6]) cancel. The unlinked diagrams must cancel in order to 
have an energy per lattice site that is finite in the thermodynamic limit. The nonvanishing 
fourth-order terms fall into three categories: those that link two distinct sites [Figure 1(b)]; 
three distinct sites [Figure 1(c)]; or four distinct sites [Figure 1(d)]. 

The matrix elements for the effective Hamiltonian are determined by introducing appro- 
priate complete sets of states between each of the operator factors in Eqs. (|J) and (0) for 
each of the possible intermediate virtual states summarized by the diagrams in Figure 1. 
In the following section, thermodynamic phase transitions are determined from the effective 
Hamiltonian via a mean-field-theory analysis. An approximation is made here, that the 
transition temperature T c is much smaller than the excitation energy of the lowest state 
above the ground state for the local anharmonic phonon with zero or two electrons on that 
site. In the harmonic case (oe a nh = 0), this excitation energy is Q, and even in the anhar- 
monic case, this excitation energy remains of the order of magnitude of Q. If the transition 
temperature is much smaller than the excitation energy, then restriction can be made to the 
degenerate subspace Q corresponding to the ground state of the phonons and one need not 
consider the effective Hamiltonian in the subspaces corresponding to excited phonon states. 
This is not a restrictive approximation, because the small-frequency limit is already known 
to be singular, since the degenerate subspace Qo becomes much larger when Q = than 
when Q ^ 0. Q/t* will be set equal to 0.5 for the numerical work in Section III. 

Determination of the effective Hamiltonian proceeds in a similar fashion to the harmonic 
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caseS'il'i. The minimum of the anharmonic potential lies at xq, x±, and £2 when there are 
zero, one, or two electrons, respectively, at a given lattice site (in the harmonic case we have 
x n = —ng/MQ 2 ). Let | + m), |m), and | — m) denote the mth anharmonic oscillator state 
centered about the origin with zero, one, and two electrons respectively, and E + (m), E(m), 
and E-(m) denote the corresponding eigenvalues under the unperturbed Hamiltonian H Q . 
The overlaps (±m|n) will need to be calculated numerically in the general case. For the 
harmonic oscillator, a simple form is found 

1 



(±m\n) = ^=e^ s (0\[a±^} m [a^T^/s/2] n \0) . (7) 
ym\n\ 

Here (a) is the harmonic-oscillator creation (annihilation) operator. In the harmonic case, 
one can use Eq. fl7|) to determine the parameters of the effective Hamiltonian analytically. 

We will concentrate, however, on the anharmonic case, where the local phonon problem 
must be solved numerically. The unperturbed Hamiltonian Hq consists of a collection of local 
phonon Hamiltonians, one for each lattice site, with a fixed number of electrons n = 0, 1, 2 
at each lattice site i. The local Hamiltonian is 

1 i 

H ioc = t^P 2 + -Mn 2 x 2 + a anh x A + gnx , (8) 

where we ignore terms that do not depend on the phonon coordinate or momentum. For 

harmonic phonons, the unit of length is (Planck's constant is set equal to one) x* : = 

Reexpressing the Hamiltonian in terms of a dimensionless distance y := x/x* yields 

O - - d 2 

H loc = -H(n,y) , H(n, y) := -— + V(n, y) , V(n, y) := y 2 + vy 4 + wny , (9) 



with v := 2a an h/ M 2 Vt 2 and w := 2gj\Jm£l? . The Schrodinger equation H(j) m = e(m)(f) m is 
then solved numerically with the Numerov algorithmic. Schrodinger's equation is cast into 
a three term recursion relation which is iterated from the far left and from the far right, and 
matched at the middle. To start the iteration requires a good initial guess for the eigenvalue. 
This is provided by the WKBJ guess ewKBj( m ) found from 

2m + 1 = - J dy^e WKB j(m) - V{y) . (10) 
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We generally solve for the thirty lowest eigenvalues and eigenvectors in the sectors with 0, 
1, or 2 electrons. 

The effective Hamiltonian is determined by evaluating all possible contributions from 
each of the nonvanishing diagrams in Figure 1. Consider first the second-order term in 
Figure 1(a). If both site % and site j are occupied by bipolarons, then the hopping matrix 
cannot connect the two sites by the Pauli principle. Similarly, if both sites are empty, there 
is no connection by the hopping term. It is only if one site is occupied by a bipolaron and 
the other site empty, that the hopping matrix can connect through a virtual state back to 
the degenerate subspace. Consider the case with a bipolaron fj, at site i and an empty site 
at site j. The hopping perturbation breaks the electron pair, with one electron hopping 
to site j and either hopping back to site % or the other electron also hops from site i to site 
j. Both processes are illustrated schematically below (with j or J, corresponding to a single 
electron at a lattice site): 

no =► t i =► no, (ii) 

Tio =► t I =► o TI , (12) 

where only one of the two possible intermediate states is shown. 

It is convenient to express the effective Hamiltonian in terms of pseudospin operators^. 
If the lattice is bipartite, so that it can be separated into A and B sublattices with nonzero 
hopping matrix elements only between sublattices A and B, then one can define pseudospin 
operators via 

Jf := (-1)^4 , -/ ; : i-ljV , ^ := l[n n + n n - 1] , (13) 

and the factor (— 1) J is 1 for the A sublattice and (—1) for the B sublattice. The pseudospin 
operators satisfy an SU(2) algebra and form a spin-^ representation in the strong-coupling 
limit. A doubly occupied site corresponds to an up spin and an empty site corresponds to 
a down spin. The matrix elements of the effective Hamiltonian (that connects site % to site 
j) satisfy 
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H 2 {i,j)\tf) = , 

H2(i,j)\ TI) = -\jf\hj)\ TI) + j)\ IT) , 
H 2 (iJ)\ IT) = IT) + TI) , 

ff 2 (»,i)|U> = o , (14) 

which is an XX Z Heisenberg antiferromagnet 

H2 = \ Y.ij^^^WJj + jt4\ +j?\ij)[Ji M J! - -]} • (is) 

Note that the summation is not restricted to i < j (the overall factor of | is introduced 
to compensate for double counting). Be careful not to confuse the pseudospin Hamiltonian 

from the original Hamiltonian that consists of bipolarons and empty sites. 

(2) (2) 

The two parameters jn and j]_ are determined by introducing complete sets of states 
into Eq. (|5|) for each of the virtual processes in Eqs. flTT| ) and (0) and employing the 
definitions given in Eq. ([14]). The results are 

.( 2)= f> 2 (-0\m)(m\-0)(+0\n)(n\ + 0) 

H ~ m %o ij 'E + (0) + E-(0) - E{m) - E{n) ' 1 1 



.( 2)= y> 2 (-0|m)(m|+0)(+0|n)(n|-0) 
J± ~Jko 13 E + (0) + E_(0) - E(m) - E(n) * 



The fourth-order terms are all evaluated in a similar fashion. They can be separated 
into three different forms H4 = H^{b) + H&(c) + H^{d) corresponding to the three different 
linked diagrams in Figure 1. The effective Hamiltonian for each of these three cases takes 
the form 



HM = lj2{j±^j)\[JtJi + Ji4\ +3\%JMJ! - \\\ ■ (is) 



1 -j 



2 . , 



+ fx^WJi + JTJf 1 + K +iHWJi - i]} , (19) 



i,j,k,l 
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h - [Ji~jj + Ji Jf + Ji J\ + J-i j? + JkJj + ^ Ji + ^fc"Jj + <4 
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+ 209 - e)[j?j;(j£jf + J k jn + J z k Jf{JtJj + JTJf) 
+ 2( 7 - ^[j? r k {jpr + JjJt) + JpiiJtJk + Ji4)\ 

-( J+ J fc " + JrJ+)(J+Jr + J/ J, + )] + [-8a + 45 + 2i/] J? J* J£ J? } , (20) 

where the lattice-site-index dependence of the parameters has been suppressed, and the 
prime on the summations means that the sites k and k, I are all distinct (the overall 
factors of \ and | are introduced to compensate for double counting). Explicit expressions 



for each of the parameters in Eqs. flTq), Q19D , and (p0|) are given in the appendix. 

The effective pseudospin Hamiltonian is an anisotropic frustrated antiferromagnetic 
Heisenberg model with additional quartic spin-spin interactions. 

The effect of the Coulomb repulsion is almost trivial. The Coulomb repulsion does not 
change any of the matrix elements, all it does is shift the energy of the bipolaron upward 
by U c : E-(m) — > E-(m) + U c , so the effect of Coulomb repulsion can be included without 
much extra effort. One might try to approximate the effect of the anharmonicity by an 
instantaneous Coulomb repulsion, chosen to match the reduction of the bipolaron binding 
energy when the anharmonicity is included, i. e. define 

U c (a anh ) :=E_(0;a anh )+E + (0;a anh )-2E(0;a anh ) + \U\ , (21) 

as the difference in the binding energy of the anharmonic bipolaron from the harmonic 
bipolaron. This reduction in bipolaron binding energy is plotted in Figures 2(a) and (b) 
as a function of the harmonic bipolaron binding energy (for fixed values of a an h) and as 
a function of a an h (for fixed values of U), respectively. Notice how even a small value of 
the anharmonicity produces a large reduction of the bipolaron binding energy, and that 
the reduction increases as the electron-phonon interaction strength increases (the phonon 
frequency is fixed at fl/t* = 0.5). 
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III. MEAN-FIELD-THEORY ANALYSIS 



We consider the case of hopping between nearest neighbors on a hypercubic lattice in 
(i-dimensions. The only nonzero matrix elements are then = t* /2\fd when i and j are 
nearest neighbors. As described above, t* is chosen to be the energy unit. This choice of 
scaling the hopping matrix elements inversely as the square root of the dimensionality is 
made so that the theory has a nontrivial limh@ as d — > oo. 

We employ a mean-field-theory analysis to the pseudospin form of the effective Hamil- 
tonian. Two types of phase transitions occur: (1) staggered order along the z-axis [corre- 
sponding to CDW order at the "antiferromagnetic" (tt, it, . . .) point]; and (2) staggered order 
along the x-axis (corresponding to SC order with a zero-momentum pair-field state). The 
mean-field theory becomes exact in high dimensions, and should provide an upper bound 
to the transition temperatures in finite dimensions (because nonlocal quantum fluctuations 
should reduce the transition temperature further). 

The mean-field theory is constructed by determining the molecular field at each lattice 
site and equating the expectation value of the magnetization with that of a free spin in an 
external magnetic field equal to the molecular field, h mo i_, yielding 

(J>= i^rjtanh^lh^.l , (22) 

I "-mol. | ^ ^ 

as first described by Gorterll. The difficult part of the calculation involves a correct determi- 
nation of the molecular field. One must be certain to properly count the contributions from 
each of the diagrams in Figure 1. Note that the total number of distinct nearest-neighbor 
pairs corresponding to Figures 1(a) and 1(b) is Nd, with N the number of lattice sites, and 
each pair appears twice in the unrestricted summations. There are two classes of second- 
neighbor diagrams corresponding to Figure 1(c): those where j and k are not parallel [of 
the (1, 1, 0, . . .) form] and those where j and k are parallel [of the (2, 0, 0, . . .) form]. In the 
first case, there are Nd(d— 1) pairs, with each pair appearing four times in the unrestricted 
summation. In the second case, there are Nd pairs, with each pair appearing twice in the 
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unrestricted summation. Finally there are ~Nd(d — 1) distinct squares corresponding to 
Figure 1(d) with each square appearing eight times in the unrestricted summation. Using 
these results, it is a straightforward exercise to rewrite the Hamiltonian as a summation 
over distinct pairs (and squares) and then extract the molecular fields for the corresponding 
ordered phases. 

A hypercubic lattice is bipartite, so it divides into two sublattices A and B, where 
the nearest-neighbor hopping occurs only from sublattice A to sublattice B or vice versa. 
The paramagnetic high-temperature phase corresponds to a uniform magnetization of the 
pseudospins on each sublattice: 

(J A ) = (J B ) =: ^me z = ^(p e - l)e z , (23) 

where e z is the unit vector along the z-axis and p e is the electron concentration. The 
self-consistent equation for the pseudospin magnetization becomes 



m 



-v) 



tanh -U)+md [-jj 2) - jj 4) + (2d - - + (d - 1)(5 + \i 

+m 3 d(d - 1) [2a - 5 - -u] | . (24) 

The dependence of the chemical potential \i upon the electron concentration p e can easily 
be determined by inverting Eq. (|24|). 

The transition temperature to the commensurate charge- density- wave phase occurs at a 
temperature where the pseudospin magnetization satisfies 

1 1 

(Ja) = -(m + m)e z , (J B ) = -(m - m )e g , (25) 

in the limit m' — > 0. We only consider the transition to a commensurate CDW, because 
transitions to incommensurate phases can only occur if the frustration induced by the fourth- 
order terms in the effective Hamiltonian becomes large enough. Since the validity of the 
truncated strong-coupling expansion fails when the fourth-order terms are comparable in 
size to the second-order terms, we ignore the complication of incommensurate order here. 
The transition temperature is then easily found to be 
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T c = \ Pe (2 - Pe ){djW - d 2 [Q 3 [ + 2j'{ -5+~v + (2a-5-~v)(l- Pe f] 

+rf[j{ 4) +3j[ | +j[ | '-5 + ^+(2a-5-^)(l-p e ) 2 ]} . (26) 

Note that this expression differs slightly (in the coefficient of the 5 term) from that given 
previously!, and corrects a typo in that work. Explicit formulae for the parameters appearing 
in Eq. fl26|) appear in the Appendix. 

Likewise, the superconducting transition temperature is determined by finding the tem- 
perature where the pseudospin magnetization satisfies 

(Jl) = (J z b) = \rn , (J%) = -<Jg) = l -m! , (27) 

in the limit m! — > 0. The transition temperature is 

Tc = i r P uo 1 m W± + ^1 ~ 2f± + 2/3 - 7 + 2e - » + (2/3 - 7 - 2e + //)(1 - Pe f] 

+d[ 3 ( l ] - 2f ± + j'i - 2/3 + 7 - 2e + /x + (-2/? + 7 + 2e - //)(1 - p e ) 2 ]} . (28) 

Explicit formulae for the parameters appearing in Eq. (pS|) appear in the Appendix. 

The above expressions for the transition temperatures to CDW or SC order are evaluated 
below in the infinite-dimensional limit. In this case, there is no contribution from the fourth- 
order terms that are multiplied by a linear power in d, because they scale like t* 4 /d — > in 
the large-dimensional limit. 

At half-filling (p e = 1), the CDW-phase is expected to be the ground state if either the 
anharmonicity or the Coulomb repulsion is not too large. Figure 3 plots the CDW transition 
temperature at half filling for different values of the anharmonicity 3(a) and the Coulomb 
repulsion 3(b). Both the second-order approximations (which monotonically diverge as the 
interaction strength approaches zero) and the fourth-order approximations (which properly 
show a peak in T c as a function of the coupling) are plotted. In Figure 3(a), four differ- 
ent values of the anharmonicity are shown: a a nh = 0.0 (the harmonic case) (solid line); 
a anh = 0.001 (dotted line); a an h = 0.005 (dashed line); and a an h = 0.015 (dash-dotted line). 
It is apparent that even though the second-order approximation shows large enhancements 
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to the transition temperature as the anharmonicity is increased, the fourth-order approxi- 
mation indicates that the maximal CDW transition temperature actually decreases as the 
anharmonicity increases. Furthermore, the value of coupling strength where the maximum 
occurs increases as a function of anharmonicity. This is to be expected since the anhar- 
monicity acts in some respects like a retarded Coulomb repulsion. What is surprising is that 
relatively small values of anharmonicity have such large effects on the transition temperature. 

In Figure 3(b), the CDW transition temperature at half filling is plotted for four different 
values of U c (a an h)- The same values of a an h are used as in Figure 3(a). Note that once again 
the maximal T c decreases as a an h increases (implying that U c increases). Furthermore, the 
peak does not move as rapidly to larger values of the coupling strength, indicating that 
the retardation effects are rather strong even at the relatively large phonon frequency of 
Q/t* = 0.5. 

Note that these strong-coupling phase diagrams, may be more accurate than those of the 
harmonic model, because it is known that in the case of either an anharmonic interaction, 
or a Coulomb repulsion, that the CDW instability only occurs when the coupling strength 
U is large enough in magnitude. Hence the true T c does vanish at a finite value of U, as it 
does in the fourth-order approximation. 

Our analysis at half filling ignored the possibility of SC order. In fact, the system will 
go superconducting at half filling in the regime where the CDW transition temperature has 
been suppressed to zero, but since this regime corresponds to a region where the fourth- 
order approximations are breaking down, we have not complicated Figure 3 by including 
the superconducting solutions at small \U\. Rather, we examine what happens as the system 
is doped away from half-filling. In Figure 4, the phase diagram is plotted for the coupling 
strength U = — 1.5625£* which lies at the peak of the CDW transition temperature curve 
when a an h = 0. This case represents a lower limit of applicability of the strong-coupling 
expansion. In Figure 4(a), the anharmonicity varies from a an h = 0, 0.003, 0.005, 0.01. 
The solid lines denote CDW solutions, and the dotted lines are SC. Note that the CDW- 
SC phase boundary lies at exponentially small densities in the harmonic case. As the 
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anharmonicity is turned on, the CDW transition temperature is suppressed, and the SC 
transition temperature is enhanced, so the CDW-SC phase boundary moves in towards half 
filling. When a an h = 0.005, the CDW phase has disappeared. The enhancement of the 
SC transition temperature continues as a an h increases, until it is significantly larger than 
the maximal CDW transition temperature of the harmonic Holstein model. Hence, the 
strong-coupling approximation, through fourth order, predicts a large enhancement in the 
SC transition temperature relative to the harmonic case. It is not clear whether this result 
is an artifact of the approximation, or is a real effect. Note further that the electron density 
where the maximal superconducting T c occurs is near the band edges. This is similar in spirit 
to Hirsch's mechanism, but, as far as we can tell is unrelated, because in the strong-coupling 
limit the effect occurs both at the top and the bottom of the band, since the approximation 
explicitly retains electron-hole symmetry. 

In Figure 4(b), the same phase diagrams are plotted, this time using U c (a an h) with the 
same values of a anh as in Figure 4(a). Clearly one can see that the effect of the Coulomb 
repulsion is quite different from the anharmonicity. The phase diagram does not change 
much, and the large enhancement of T c in the SC channel does not occur. 

In order to check to see whether these results are generic, or occur simply because one 
is close to the limiting region where the approximations are expected to hold, we have also 
calculated the phase diagrams for a stronger value of the interaction strength U = —AM*. 
In Figure 5(a), four different values of a anh (0.0, 0.01, 0.03, 0.05) are plotted. Here the 
results are similar to those found in Figure 4, except the anharmonicity initially causes the 
CDW T c to rise, because of the reduction in the bipolaron binding energy. The SC transition 
temperature still has a large maximum at low densities, and is significantly enhanced relative 
to the harmonic case. Furthermore, the CDW-SC phase boundary continues to move toward 
half-filling until it disappears at a an h ~ 0.03. In Figure 5(b), we plot the same phase 
diagrams, this time using U c (a a nh)- Once again, the phase diagram displays very different 
behavior, with the CDW-SC phase transition remaining at exponentially small densities. 

Finally, we study how the critical electron density, where the CDW-SC phase boundary 
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lies, evolves as functions of U and a an h- In the anharmonic Holstein model, we see in 
Figure 6(a) that the phase boundary moves very rapidly as the anharmonicity is turned 
on. This indicates how the anharmonicity strongly favors SC solutions relative to CDW 
order. In Figure 6(b), we show the analogous plots of the critical electron density for the 
harmonic Holstein-Hubbard model, with U c chosen from Eq. QZTD, and the same values of 
otanh- Clearly the anharmonic behavior is not easily mimicked by an instantaneous Coulomb 
repulsion, and the retardation effects cannot be neglected. 

IV. CONCLUSIONS 

The strong-coupling expansion for the anharmonic Holstein-Hubbard model has been 
presented through fourth order in the hopping. This result extends the analysis of the har- 
monic easel. We find some interesting results from this analysis. First, the anharmonicity 
reduces the bipolaron binding energy and second, it reduces the the equilibrium phonon- 
coordinate spacing between one electron at a site and either zero or two electrons at a 
site. The first effect is expected to cause an enhancement to transition temperatures in the 
strong-coupling regime, and the second should enhance the Franck- Condon overlap factors 
for superconducting order, favoring the SC phase relative to the CDW. Similarly, a Coulomb 
repulsion will reduce the bipolaron binding energy, but does not alter the Franck-Condon 
overlaps. We find that in the CDW phase, the maximal transition temperature actually 
decreases when either anharmonicity or Coulomb repulsion are turned on. For the super- 
conducting phase, the enhancement of the Franck-Condon overlap factors (equivalent to a 
widening of the polaron band) causes a large enhancement of the SC transition temperature 
at low electron density, even for moderate values of the anharmonicity. In no case do we find 
that the effect of the anharmonicity is easily mimicked by an effective Coulomb repulsion. 

Since the strong-coupling expansion is expected to fail in both the low-electron- 
concentration regime, and when the effective bipolaron binding energy is no longer much 
larger than the hopping integral, it is possible that this large enhancement of the super- 
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conducting transition temperature is just an artifact of the current approximation. It is 
important to compare these approximations to exact quantum Monte Carlo simulations of 
the transition temperatures of the anharmonic Holstein model. Work in this direction is 
currently in progress. 
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APPENDIX A: APPENDIX: PARAMETERS OF THE FOURTH-ORDER 
EFFECTIVE PSEUDOSPIN HAMILTONIAN 

In this appendix explicit expressions are given for the parameters that appear in the 
fourth-order effective pseudospin Hamiltonian, summarized in Eqs. (H), fll9f) , and (p0|) . 
The notation used is that given in the text. 

First the parameters in Eq. ([18|) : 



(+0\n)(n\ +m)(+m\l)(l\ + 0)(-0\n')(n'\ - m'){-m'\l')(l'\ - 0) 
. [E+(0) + E_(0) - E{n) - E(n')][E+(0) + E_(0) - E+(m) - E_(m')] 



E+(0) + £L(0) -E(l) -E{V) 
(+0\n)(n\ - m)(-m\l)(l\ + 0}(-0|n')(n'| + m')(+m'\l')(l'\ - 0) 

[E+(0) + £_(0) - E(n) - E(n')}[E + (0) + EjO) - E+(m) - E_(m')} 



ACKNOWLEDGMENTS 




oo 



J 



771,777/ ,n, n f —0 
m-(- m' ^0 



1 



X 



E+{0) + £L(0) -E(l) -E{V) 
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- E 



{+Q\n){n\ + 0)(-0|n')<n'| - 0)(+0|/)(/| + 0)(-0|/')(/'| - 0) 



,^=oU^+(°) + ^-(°) " E ( n ) ~ E(n')) 2 [E + (0) + E.(0) - E(l) - E(l>)} 
(+0\n)(n\ - 0)(-0|n')(n / | + 0)(-0|/)(/| + 0)(+0|/')(/'| - 0) 



+ 



[£+(0) + £L(0) - E(n) - E(n')} 2 [E + (0) + £_(0) - E(l) - E(l')} 



(Al) 



A 4) 



-8t] 



E 



l ,V ,rn,m f ,n,n f =0 



+ 



[+0\n)(n\+m)(+m\l)(l\ - 0) (-0|n') (n'\ - m')(-m'\l')(l'\ +0) 

L [E + (0) + E_(0) - E{n) - E(n')][E+(0) + £_(0) - E+(m) - Ejw/j] 

1 

* E+(0) + E_(0) - E(l) - E(l') 
(+Q\n)(n\ -m)(-m\l)(l\ - 0)(-0|n')(n'| + m'){+m'\l'){l'\ +0) 

[E+(0) + E_(0) - E{n) - E(n')][E+(0) + E_(0) - E + {m) - E_(m')} 

1 

X E+{0)+E_{0)-E{l)-E{l'). 

{+0\n)(n\ + 0)(-0|n / )(w / | - 0)(+0|/)(/| - 0)(-0\l')(l'\ + 0) 
[E + (Q) + £_(0) - E(n) - E(n')} 2 [E + (0) + E.(0) - E(l) - E(l')} 
(+0|n)(n| - 0)(-0\n')(n'\ + 0)(-0|/)(/| - 0) (+0|Q (/'| + 0) 
[E+(0) + £L(0) - £(n) - E(n')} 2 [E + (0) + £_(0) - £7(0 - £(/')] 

Next the parameters in Eq. (|T 



- E 

l,l',n,n'=Q 
+ 



(A2) 



4l 



E 

Z,Z',m,n,n'=0 

(+0|n)(n| + 0)(-0|n / )(n'| + m)(+m\l')(l'\ - 0)(+0|/)(/| + 0) 

[E+(0) + £_(0) - E(n) - E(n')][2E+(0) + £_(0) - £+(m) - £(ra) - E(l)] 

1 1 



x 



E + (Q) + E_(Q) -E(l) -E{1') E+(0) + E-(0) - E(n) - E(V) 
\+0\n)(n\ + 0)(-0\n')(n'\ - m)(-m\l')(l'\ - 0)(+0|/)(/| + 0) 

u'^o l E +(°) + ^-(0) - E(n) - E(n')][E + (0) - E + (m)][E + (0) + E_(0) - E(l) - E(l')] 
(+0\n)(n\+0)(-0\n')(n'\-0)(-0\l')(l'\-0)(+0\l)(l\ + 0) 



+2 E 



mjtQ 
OO 



\,'£=o [E+(0) + J5_(0) - £(n) - £(n')] 2 [£+(0) + J5_(0) - E(l) - E{V)\ 



(A3) 



j'{ - ~^ t % t % E 

1,1' ,m,n,n'=0 



X 



(+0|n)(ra| + m)(+m|/ / )(/ / | + 0) (-0\n') (n' '\ - 0)(+0|/)(l| + 0) 
E + (0) + £7_(0) - E(n) - E(n') 
1 

[2£7+(0) + £_(0) - E+(m) - E(l) - E(n')][E+(0) + £L(0) - E(V) - E(n')} 



(A4) 



, OO 



1,1' ,m,n,n'=0 
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(+0\n) (n\ + 0) (-0\n') (n'\ + m) (+m\l') + 0) (+0|Z) (/| - 0) 

[£+(0) +E-(0) -E(n) - E(n')][2E+(0) + E_(0) - E + {m) - E{n) - E(l)] 

1 1 

X — =777T + 



+2 E 



111 



E 



E+(0) +E-(0) -E(l) -E(l') E+(0) + E-(0) - E(n) - E(l') 
(+0\n)(n\ + 0)(-0\n')(n'\ - m){-m\l')(l'\ + 0)(+0|Z)(Z|- 0) 

tX'=o l E +(°) + E -(°) - W) - E(n')}[E + (0) - E + (m)}[E + (0) + KjO) - E(l) - E{1')} 
(+0\n)(n\ + 0)(-0|n')(n'| - 0)(-0|/')(/'| + 0)(+0|/)(/| - 0) 



m^O 
OO 



M ,^, =0 [JS? + (0) + J5_(0) - E(n) - E(n')][E + (0) + E.(0) - E(l) - E(l')} 

1 1 



x 



+ 



E+(0) + E_(0) - E(n) - E(n') E+(0) + £L(0) - E{1) - E{V) 



(A5) 



{OO 



+2 E 



i,i',m,n,n'=0 

(+0|n)(n| + m)(+m|Z')(Z'| + 0}(-0|n')(n / | + 0)(+0|Z)(Z| - 0) 

[E+(0) + E-(Q) - E{n) - E(n')][2E+(0) + £_(0) - E+(m) - E(ri) - E(l)} 

1 

X E+{0) + E_(Q) - E(l) - E{V) 
(+0|ra)(w| - m)(-m\r)(l'\ + 0)(-0\n'){n'\ + 0)(+0|Z)(Z| - 0) 

,^ ;=0 [E + (0) + £_(0) - E(n) - E(n')][E + (0) - E + (m)][E + (0) + E_(0) - E{1) - E{1')\ 
(+0\n)(n\ - 0)(-0\n')(n'\ + 0)(-0\l')(l'\ + 0)(+0|/)(/| - 0) 



oo 



a 



2 M ,5 =0 [^+(0) + E-(0) - E(n) - E(n')nE + (0) + E-(0) - E{1) - £(Z')] 
Finally the parameters in Eq. (p]|): 

(+0|n)(n| + 0)(-0|n')(n'| - 0>(+0|Z)(Z| + 0)(+0|Z')(Z'| + 0) 



^tijtjktkltli ^ ] 



M ,^, =0 [£+(0) + tf_(0) - £(n) - E(n')][E + (0) + tf_(0) - £(Z) - E{n')\ 
X £+(0) + Ejfi) - EiV) - E(n>) ' (A7) 



(A6) 



(+0|n)(n| - 0}(-0|n')(n / | + 0>(+0|Z)(Z| + 0>(+0|Z')(Z'| + 0) 



/3 2 W»*« ^g^l [£+(0) + E-(0) - E(n) - E(n>))[E + (0) + E.(0) - E(l) - E{n)\ 



x- 



x 



£+(0)+EL(0) -E{V) -E{n) 
(+0\n){n\ + 0)(-0|n')(n'| + 0)(+0|Z)(Z| + 0)(+0|Z')(Z'| - 0) 

E+(0) + E_(0) - E(n) - E(n') 
1 

[E + (0) + E.(0) - W) ~ E(n)][E + (0) + £L(0) - £(Z) - E{V)} 

1 

+ [£+(0) + EJ0) - E(l) - E(n')}[E + (0) + E_(0) - £(Z) - E(l')} 

1 



+ 



[£+(0) + £_(0) - £(Z) - £(n')p+(0) + £_(0) - E(Z') - E{n>) 
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« f (+0\n)(n\ + 0)(-0|n')(n'| + 0)(+0|Z)(Z| - 0)(+0\l')(l'\ + 0) 
7 4 ^W« 2.^| [E+(Q) + E _ (Q) _ E{n) _ E (n')][E + (0) + E40) - E(l) - E(n>) 



1 



x- 



+ 



X 



E + (0) + E_(0) - E(l) - E(l>) 
(+0\n)(n\ + 0)(-0\n')(n'\ + 0)(+0|Q(/| + 0)(+0\l')(l'\ - 0) 
[E+(0) + E_(0) - E(n) - E(n')][E+(0) + E_(0) - E(l) - E(n)} 

1 1 

+ 



5 = 



E+(0) + £_(0) -E(V) -E{n) E+(0) + £_(0) - E(l) - E(l') 

~ f (+0|n)(n| + 0)(-0K)(n-|-0)(+0|/)(/| + 0)(-0|Q(n-0) 
- % i^^t [E+(0) + E _ (0) _ E{n) _ E (n>)][E + (0) + E.(0) - E(l) - E(n>) 



, (A9) 



1 



x- 



+ 



E+(0) + £_(0) -E(l) -E{V) 
(+0\n){n\ + 0)(-0|n')(n'| - Q> <— Q|Z> <Z| - 0)(+0|/ / )(/ / | + 0) 

[E+(0) + £_(0) - £(n) - E(n')][E+(0) + £_(0) - E(l) - E(n)\ 

1 

x 



e — 2tijtjktkitu ^2 I 2 



x 



£+(0)+£_(0)-£(Z) -£(Z'). 

(+0|n)(n| - 0)(-0\n')(n'\ + 0)(+0|Z)(Z| + 0) (-0|Z') (Z'| - 0) 

E+(0) + f?_(0) - £(ra) - f?(n' 
1 

X 2£+(0) + 2£_(0) - E{n) - E(n') - E(l) - E{1') 

1 1 

+ 



+ 



E + (0) + E-(0) - E(l') - E(n') E+(0) + £_(0) - E(l) - E(n) 
(+0\n)(n\ + 0)(-0|n / )(n / | - 0)(+0|Z)(Z| - 0)(-0|O(Z'| + 0) 

[E+(0) + E-(0) - E(n) - E(n')][E+(0) + £_(0) - E(l) - E(n')\ 
1 1 



x 



+ 



+ 



E+(0) + E_(0) - E(l') - E(n') E+(0) + E_(0) - E(l) - E(l>) 
+%)H + 0)(-0|n')(n'| - 0)(-0|/)(/| + 0)(+0|Z')(Z'| - 0) 

[E+(0) + EJO) - E(n) - E(n')][E+(0) + E_(0) - E(l) - E(n 

1 1 

+ 



x 



A* — —2tijtj k tkitii ^ 1 ^ 



E+(0) + £_(0) - £(Z') - E(n) E+(0) + £_(0) - E(l) - E(V) _ 
(+0\n)(n\ - 0)(-0\n')(n'\ - 0)(+0|/)(/| + 0)(-0|Z')(/'| + 0) 



E+(0) + £_(0) - E(n) - E(n' 
1 



x 



2E + {0) + 2£_(0) - E{n) - E(n') - E(l) - E(V) 



X 


1 




1 


E + {0) + E_(0) - E{1>) - 


£(n') 


' £ + (0) + £_(0) -£(Z) 




(+0|n)(n| -0)(-0|n')(n'| 


+ 0)(- 


-0|Z)(Z| + 0)(+0|Z / )(Z , | + 0) 



[£+(0) + £_(0) - £(ra) - £(n')P+( ) + E_(0) - f?(Z) - E(n)\ 
1 

x- 



E+(0) + E_(0) - E(Z') - 

23 



(A10) 



, (All) 



(+0|n)(n| + 0)(-0|w')(n'| + 0)(+0|Z)(Z| - 0)(-0|Z')(Z'| - 0) 

+ [E+(0) + £_(0) - E(n) - E(n')][E+(0) + £_(0) - E(l) - E(n')\ 
1 

X E+(0) + EJO) - E(l') - E{n') 



(A12) 



v — —stijtjktkitu 



(+0\n){n\ + 0)(-0\n')(n'\ - 0)(+0|Z)(Z| + 0)(-0\l')(l'\ - 0) 
E+(0) + E_(0) - E(n) - E(n' 
1 



x 



x 



2E + (0) + 2E_(0) - E(n) - E(n') - E{1) - E(V) 

1 1 

+ 



£+(0) +£_(0) -E(l') -E(n) E+(0) + £_(0) - E(l) - E(n') 



, (A13) 



P = 



u , , , ^ f (+Q|») ("I - Q) (-%') fo'l + Q) (+Q10 (-Q10 + Q) 



£+(0) + £_(0) - E(n) - E(n' 
1 



x 



X 



2E + (0) + 2E_(0) - E{n) - E{n') - E(l) - E(V) 
1 1 



+ 



+ 



E+{0) + E_(0) - E{1) - E{V) E + (0) + E_(0) - E(n) - E(n') 
(+0|n)(n| - 0)(-0\n')(n'\ + 0)(-0|/)(/| + 0)(+0|Z')(Z'| - 0) 

[E + {0) + Ejfi) - E(n) - E(n')][E + (0) + EjO) - E{1) - E{n')] 

1 1 

+ 



x 



+ 



E+(0) + E-(0) - E(l) - E{1>) E + (0) + £_(0) - E(l') - E(n') 
(+0\n}{n\ - 0)(-0|w')(w'| + 0)(+0|Z)(Z| - 0)(-0|Z')(Z'| + 0) 

[E+(0) + E-(0) -E(n) - E(n')][E+(0) + E-(0) - E(l) - E(n)\ 
1 1 

X -77^77 77-777 77777 777"^ + 



E+(0) + E_(0) - E{V) - E(n) E+(0) + E_(0) - E(l) - E(V) 



• (A14) 
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FIGURES 

FIG. 1. Schematic diagrams used in the determination of the effective Hamiltonian. The 
second-order diagram is plotted in (a), while the fourth-order diagrams appear in (b-e). The 
fourth-order diagrams link two (b), three (c), or four (d) distinct lattice sites. The contributions 
from the unlinked diagram (e) vanish. 

FIG. 2. Reduction of the bipolaron binding energy [which is the definition of U c (a an h)] plotted 
as functions of the electron-phonon coupling (a) and the anharmonicity (b). In Fig. 2(a), three 
different values of the anharmonicity parameter are chosen: a an h = 0.001 (solid line); 0.01 (dotted 
line); and 0.1 (dashed line). Note how the reduction of the bipolaron binding energy increases with 
both a an h and U. In Fig. 2(b), the reduction in bipolaron binding energy is plotted as a function 
of a an h for three values of \U\: \U\ = 1 (solid line); \U\ = 4 (dotted line); and \U\ = 16 (dashed 
line). Note that even relatively small values of a an h produce a sharp reduction in the bipolaron 
binding energy. 

FIG. 3. Charge-density-wave transition temperature at half filling for the anharmonic Holstein 
model (a) or the Holstein-Hubbard model (b). Both second-order and fourth-order approximations 
are plotted (the fourth-order approximations have a peak in T c ). Note that the maximal T c for 
the fourth-order calculation decreases as either the anharmonicity or Coulomb repulsion is turned 
on. In Fig. 3(a), the anharmonicity varies from a an h = 0.0 (solid line), a an h = 0.001 (dotted 
line), a an h = 0.005 (dashed line), and a an h = 0.015 (dash-dotted line). In Fig. 3(b), the Coulomb 
repulsion is chosen to match the reduction in the bipolaron binding energy for each value of the 
anharmonicity plotted in (a), i. e., U c = Ucis^anh) with a an h = 0.0, 0.001, 0.005, 0.015. Note that 
the effect of anharmonicity is not easily mimicked by Coulomb repulsion except for the smallest 
values of a anh . 



27 



FIG. 4. Phase diagram of the anharmonic Holstein model (a) and the Holstein-Hubbard model 
(b) for \U\ = 1.5625 corresponding to the peak in the charge-density-wave transition tempera- 
ture at half filling for the harmonic Holstein model. The solid lines denote charge-density-wave 
transition temperatures and the dotted lines denote superconducting transition temperatures. In 
Fig. 4(a), the anharmonicity strength assumes four values ctanh — 

0.0, 0.003, 0.005, 0.01. Note that 

the anharmonicity strongly enhances superconducting solutions relative to charge-density-wave so- 
lutions, and if a an h is large enough, there is no CDW order. Note furthermore, that the maximal 
superconducting T c is larger than the maximal CDW T c in the harmonic case and occurs at a 
rather low value of the electron filling. In Fig. 4(b), the Coulomb repulsion is fixed to be U c (a an h) 
with a an h the same as in (a). Note that the effect of the Coulomb repulsion on the transition 
temperatures is relatively minor, and the strong enhancement of the superconducting T c does not 
occur. 

FIG. 5. Phase diagram of the anharmonic Holstein model (a) and the Holstein-Hubbard model 
(b) for \U\ = 4.0. The solid lines denote charge-density- wave transition temperatures and the 
dotted lines denote superconducting transition temperatures. In Fig. 5(a), the anharmonicity 
strength assumes four values a an h = 0.0, 0.01, 0.03, 0.05. Note that the large enhancement in the 
superconducting T c is seen here too (although it occurs at a larger value of the anharmonicity). In 
Fig. 5(b), the Coulomb repulsion is fixed to be U c (a an h) with a an h the same as in (a). Once again, 
the effect of Coulomb repulsion is drastically different from that of the anharmonicity. 

FIG. 6. Plot of the charge-density-wave-superconductor phase boundary as functions of elec- 
tron concentration and U. In Fig. 6(a), the phase diagram is plotted for the anharmonic Holstein 
model with a an h = 0.0 (solid line), 0.001 (dotted line), 0.005 (dashed line), and 0.015 (dash-dotted 
line). In Fig. 6(b), the phase diagram is plotted for the Holstein-Hubbard model with U c = U c (a an h) 
and a an h chosen to have the same values as in (a). These figures show that the anharmonicity 
strongly favors superconductivity relative to CDW order, much more so than turning on a Coulomb 
repulsion. 
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